
#include "CFDPassFlowIC.h"
#include "CFDProblem.h"
#include "Attitude.h"
// #include "Eigen/Geometry"
// using namespace Eigen;

template<>
InputParameters validParams<CFDPassFlowIC>()
{
  InputParameters params = validParams<CFDInitialCondition>();
  params.addParam<Real>("velocity", 1, "初始的速度");
  return params;
}

CFDPassFlowIC::CFDPassFlowIC(const InputParameters & parameters) :
	CFDInitialCondition(parameters),
	_velocity(getParam<Real>("velocity")),
	_attitude(_cfd_problem._attitude),
	_gamma(_cfd_problem._gamma),
	_mach(_cfd_problem._mach)
{}

void CFDPassFlowIC::compute()
{
	InitialCondition::compute();

	if(!_constant_ic)
		return;

	NumericVector<Number> & solution = _var.sys().solution();
	std::vector<dof_id_type> dof_indices;
	dof_indices = _var.dofIndices();

	Number fineval = value(_current_elem->centroid());
	solution.set(dof_indices[0], fineval);
	const unsigned int n_dofs = dof_indices.size();
	for (unsigned int i = 1; i < n_dofs; i++)
	{
		solution.set(dof_indices[i], 0);
		_var.setNodalValue(0, i);
	}
}

Real CFDPassFlowIC::value(const Point & p)
{
	Real density = 1;
	Vector3d vel  = _velocity*(_attitude.earthFromWind()*Vector3d::UnitX());
	if(_current_elem->dim() < 3)
		vel(2) = 0;
	if(_current_elem->dim() < 2)
		vel(1) = 0;
	Vector3d mom = density*vel;

	Real pre = 1./_gamma/_mach/_mach;
	Real rhoe = pre/(_gamma-1) + 0.5*density*(_velocity*_velocity);
	_cfd_data.uh[0] = density;
	_cfd_data.uh[1] = mom(0);
	_cfd_data.uh[2] = mom(1);
	_cfd_data.uh[3] = mom(2);
	_cfd_data.uh[4] = rhoe;

	return _cfd_data.uh[_component];
}
